rm(list = ls())
library(sf)        
library(haven)
library(corrplot)
library(osmdata)    
library(ggplot2)    
library(tidyverse)
library(fixest)
library(sandwich)
library(DirectEffects)
# install.packages("remotes")
# remotes::install_github("grantmcdermott/ritest")
library(ritest)


# Load Pre-processed Data --------------------------------------------

df <- read_dta("data/clean_survey.dta")

# Figure 1 ----------------------------------------------------------

# load and project satellite detected damage 
damage <- st_read("data/raw_data/UNITAR/Mosul_Damage_Sites_20170804.shp")
st_crs(damage$geometry) <- 4326
st_transform(damage$geometry, crs = 4326)

# using the osmdata package, load and project streets and rivers in Mosul from OpenStreetMap 
Mosul <- getbb("Mosul Iraq")
Mosul[1,] <- c(43.05, 43.24)
Mosul[2,] <- c(36.29, 36.41)
streets <- Mosul%>%
  opq()%>%
  add_osm_feature(key = "highway", 
                  value = c("motorway", "primary", 
                            "secondary", "tertiary")) %>%
  osmdata_sf()
st_crs(streets$osm_lines) <- 4326
small_streets <- Mosul%>%
  opq()%>%
  add_osm_feature(key = "highway", 
                  value = c("residential", "living_street",
                            "unclassified",
                            "service", "footway")) %>%
  osmdata_sf()
st_crs(small_streets$osm_lines) <- 4326
river <- getbb("Mosul Iraq")%>%
  opq()%>%
  add_osm_feature(key = "waterway", value = "river") %>%
  osmdata_sf()
st_crs(river$osm_lines) <- 4326
river$osm_lines <- river$osm_lines[river$osm_lines$`name:hsb` %in% "Tigris",]

# plot
ggplot() +
  geom_sf(data = streets$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .4,
          alpha = .6) +
  geom_sf(data = small_streets$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .2,
          alpha = .4) +
  geom_sf(data = river$osm_lines,
          inherit.aes = FALSE,
          color = "steelblue",
          size = 1,
          alpha = .8)  +
  geom_sf(data = damage$geometry,
          inherit.aes = FALSE,
          color = "darkred",
          shape = 4,
          size = .8,
          alpha = 1) +
  coord_sf(xlim = c(43.05, 43.24),
           ylim = c(36.29, 36.41),
           expand = FALSE) +
  theme_void()

# Figure 3 ---------------------------------------------------------

# load in anonymized respondents’ sampling coordinates with random error terms
survey_points <- read_csv("data/raw_data/survey_points.csv")

# project 
sf <- st_as_sf(survey_points, coords = c("Longitude", "Latitude"), crs = 4326, agr = "constant")

# plot survey points  
ggplot() +
  geom_sf(data = streets$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .4,
          alpha = .6) +
  geom_sf(data = small_streets$osm_lines,
          inherit.aes = FALSE,
          color = "black",
          size = .2,
          alpha = .4) +
  geom_sf(data = river$osm_lines,
          inherit.aes = FALSE,
          color = "steelblue",
          size = 5,
          alpha = 1)  +
  geom_sf(data = sf$geometry,
          inherit.aes = FALSE,
          color = "green",
          shape = 16, 
          size = .6,
          alpha = .5) +
  coord_sf(xlim = c(43.05, 43.24),
           ylim = c(36.29, 36.41),
           expand = FALSE) +
  theme_void()


# Figure A-3 ----------------------------------------------------------

# Run sequential g; Table 4 (Column 2)
m <- military_legitimacy ~ 
  # first stage 
  treated + factor(education) + age +  vote + pop_density + road_density + res_density + 
  IS_any_harm + factor(income_pre_IS) |
  # intermediate confounders 
  factor(identity) + factor(sharia) + factor(prayer) + daesh_govern + factor(daesh_corrupt) + factor(daesh_taxes)  + factor(corrupt_current) |
  # de-mediated variables 
  death_or_injury
direct <- sequential_g(m, data = df)

# plot sensitivity
out_sens <- cdesens(direct, var = "treated")
plot(out_sens)

# Figure A-4 ----------------------------------------------------------
# load health data 
data <- read_dta("data/processed_data/health_survey.dta")

# aggregate to east and west 
east_sum <- data %>%
  select(side, death, death_airstrike, death_explosions, death_gunshot, death_carbomb, death_other_conflict) %>% 
  filter(side == "East/left") %>% 
  select(-"side") %>%
  summarise_all(sum)
west_sum <- data %>%
  select(side, death, death_airstrike, death_explosions, death_gunshot, death_carbomb, death_other_conflict) %>% 
  filter(side == "West/right") %>% 
  select(-"side") %>%
  summarise_all(sum)

# summarize by proportion  
death_summary <- data.frame(
  side = factor(c("West","West","West","West","West", "East","East","East","East","East"), levels=c("West","East")),
  cause = factor( c("Airstrike","Explosion","Gunshot","Carbomb", "Other","Airstrike","Explosion","Gunshot","Carbomb", "Other"), levels=c("Airstrike","Explosion","Gunshot","Carbomb", "Other")), 
  percentage = c(west_sum$death_airstrike/ west_sum$death, west_sum$death_explosions/ west_sum$death, west_sum$death_gunshot/ west_sum$death,west_sum$death_carbomb/ west_sum$death,west_sum$death_other_conflict/ west_sum$death,
                 east_sum$death_airstrike/ east_sum$death, east_sum$death_explosions/ east_sum$death, east_sum$death_gunshot/ east_sum$death,east_sum$death_carbomb/ east_sum$death,east_sum$death_other_conflict/ east_sum$death)
)

# plot 
death_plot <- ggplot(data=death_summary, aes(x=cause, y=percentage, fill=side)) +
  geom_bar(stat="identity", position=position_dodge(), colour="black") +
  xlab("Cause of Death") + ylab("Ratio of Total Deaths")


# Figure A-5 ----------------------------------------------------------

df_cor <- df %>% 
  select(us_legitimacy, cts_legitimacy,army_legitimacy,police_legitimacy, pmf_legitimacy)%>% 
  drop_na()
cor_matrix <- cor(df_cor)

# Define custom labels for  variables
custom_labels <- c("US", "CTS", "Army", "Police", "PMF")

# Assign custom labels to row and column names of the correlation matrix
rownames(cor_matrix) <- custom_labels
colnames(cor_matrix) <- custom_labels
corrplot(cor_matrix, method = 'number', type = 'upper',tl.srt = 45, tl.col = "black") # colorful number

# Figure A-6 ----------------------------------------------------------

# run main model (Table 3 Column 4) 
base <- lm(military_legitimacy ~ treated + gender + age + factor(income_pre_IS) + 
          factor(identity) + factor(education) + vote + gov_grievance_police + 
          gov_grievance_arrested + gov_grievance_sunni + gov_grievance_protest + 
          IS_any_harm + IS_electric + IS_water + IS_zakat + pop_density + 
          road_density + res_density, data = df)

# run model with conley standard errors (requires respondent coordinates)
# conley <- feols(military_legitimacy ~  treated + gender + age + factor(income_pre_IS) +
#                   factor(identity) + factor(education) + vote + gov_grievance_police +
#                   gov_grievance_arrested + gov_grievance_sunni + gov_grievance_protest +
#                   IS_any_harm + IS_electric + IS_water + IS_zakat + pop_density +
#                   road_density + res_density, data = df, vcov = "conley")
# extract se employing various approaches  
clusterse <- list(sqrt(diag(vcovCL(base, cluster = ~ Neighborhood, type = "HC1"))))
clusterse_treat <-  list(sqrt(diag(vcovCL(base, cluster = ~ treated, type = "HC1"))))
# conleyse <- list(sqrt(diag(vcov(conley))))
bootsrapse <-  list(sqrt(diag(vcovBS(base, cluster = ~ Neighborhood))))
bootsraps_treat <-  list(sqrt(diag(vcovBS(base, cluster = ~ treated))))

# Extract coefficients and standard errors
coeff_data <- data.frame(
  calculation = c("No Adjustment", "HC1 Robust (Neighborhood Clustered)", "HC1 Robust (Side of River Clustered)", 
                  "Bootsrap (Neighborhood Clustered)", "Bootsrap (Side of River Clustered)"
                  # ,"Conley"
                  ),
  estimate = rep(summary(base)$coefficients["treated","Estimate"], 5), 
  std_error = c(summary(base)$coefficients["treated","Std. Error"],
                clusterse[[1]]["treated"], 
                clusterse_treat[[1]]["treated"], 
                bootsrapse[[1]]["treated"], 
                bootsraps_treat[[1]]["treated"] 
                # ,conleyse[[1]]["treated"]
                )
)
coeff_data$conf.low <- coeff_data$estimate - 1.96 * coeff_data$std_error 
coeff_data$conf.high <- coeff_data$estimate + 1.96 * coeff_data$std_error 

# Plot 
se <- ggplot(data = coeff_data, aes(estimate, calculation)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
  # coord_cartesian(xlim = c(-.1, .5)) +
  theme_bw() +
  theme(strip.background = element_blank(),
        axis.title.y = element_blank()) + 
  labs(x = "Estimate",
       y = "SE Calucation") 

# Figure A-7 ----------------------------------------------------------
base <- feols(military_legitimacy  ~ treated +
              gender + age + income_pre_IS_2 +income_pre_IS_3 + income_pre_IS_4 + 
               identity_2 + identity_3 + identity_4 + identity_5 + identity_6 + identity_7 + 
               education_2 + education_3 + education_4 + education_5 + education_6 + education_7 + 
               vote + gov_grievance_police +
               gov_grievance_arrested + gov_grievance_sunni + gov_grievance_protest + 
               IS_any_harm + IS_electric + IS_water + IS_zakat + pop_density + 
               road_density + res_density, data = df)

est_ri = ritest(base, 'treated', cluster='Neighborhood', level = 0.95,
                reps=10000 , seed=1234)

ggplot(data.frame(betas = est_ri$betas), aes(betas)) + 
  geom_density(fill = "blue", alpha = 0.5) +
  geom_vline(aes(xintercept = base$coefficients[2]), color = "red", linetype = "dashed", size = 1) +
  theme_minimal()+
  xlim(-.3,.3) +
  labs(x = "Treatment Point Estimate",
       y = "Density") 
table(base$coefficients[2] > abs(est_ri$betas))
table(base$coefficients[2] > est_ri$betas)

est_ri_no_cluster = ritest(base, 'treated', level = 0.95,
                           reps=10000 , seed=1234)
ggplot(data.frame(betas = est_ri_no_cluster$betas), aes(betas)) + 
  geom_density(fill = "blue", alpha = 0.5) +
  geom_vline(aes(xintercept = base$coefficients[2]), color = "red", linetype = "dashed", size = 1) +
  theme_minimal()+
  xlim(-.3,.3) +
  labs(x = "Treatment Point Estimate",
       y = "Density") 




